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In many real-world oscillator systems, the phase response curves are highly heterogeneous. How- 
ever, dynamics of heterogeneous oscillator networks has not been seriously addressed. We propose a 
theoretical framework to analyze such a system by dealing explicitly with the heterogeneous phase 
response curves. We develop a novel method to solve the self-consistent equations for order param- 
eters by using formal complex- valued phase variables, and apply our theory to networks of in vitro 
cortical neurons. We find a novel state transition that is not observed in previous oscillator network 
models. 
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Synchronization phenomena are ubiquitous in nonlin- 
ear dynamical systems, such as Josephson junction arrays 
[l[ , laser arrays and biological systems 0, 0, Q . Syn- 
chronous firing of cortical neurons is considered to play 
an active role in cognitive functions [3|, and is governed 
by the intrinsic properties of neurons as well as by the 
network connectivity. These properties include the phase 
response curve (PRC) 0, [1, @] > which describes how the 
timingof a succeeding output spike is shifted by an input 

spike 0, mm, (Mini- 

In general, we can categorize the 
phase responses of cortical neurons into two types. Type- 
I PRC has only positive values (corresponding to phase 
advances), while type II has both positive and negative 
values (corresponding to phase delays) depending on the 
phase at which a stimulus is applied [8| . Mutual synchro- 
nization of excitatory neurons may be easier with type-II 
PRC than with type-I PRC, if the PRCs of the neurons 
are homogeneous [8|. However, the PRCs recorded from 
various brain areas, which include the hippocampus 1121] . 
the entorhinal cortex lla |. t he somatosensory cortex [14| 
and the motor cortex [3, [l5| , have revealed that the PRC 
type of pyramidal neurons is highly heterogeneous, espe- 
cially if they belong to different cortical layers [ll| . Even 
if two neurons have the same PRC type, the shape of 
PRC varies significantly from neuron to neuron (see Fig. 

In general, the heterogeneity of PRCs disturbs the sta- 
bility of the synchronous state. However, the hetero- 
geneity of PRCs can be compensated by other intrinsic 
properties that enhance synchronization. In this paper, 
we explore how such compensation may occur in net- 
works of heterogeneous oscillators. The population dy- 
namics of the heterogeneous phase oscillators was first 
studied in the Kuramoto model @, which demonstrated 
the emergence of transitions between synchronized and 
desynchronized states [ljl [13, LH, Ell- These studies 
transformed the heterogeneity of the PRC shapes into 
that of the natural frequencies, assuming that the het- 
erogeneity is weak in both cases. Here, we develop an 



analytical method to explicitly deal with the heterogene- 
ity since the PRC shapes of cortical neural oscillators are 
strongly heterogeneous. Unlike in the original Kuramoto 
model, we can show that the order parameter of synchro- 
nization changes discontinuously in the neural population 
with heterogeneous PRCs. This implies that this type of 
heterogeneity creates a dramatic effect on networks of 
neural oscillators. 
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FIG. 1: Heterogeneity of the phase response curves (PRCs) 
of cortical neurons in our in vitro recording studies [lg]- We 
performed the least-square-error fitting of the PRCs with 
Z(9) = — cos((9 — an) + cos air by changing the shape param- 
eter a (see also Eq. (|lip ). (a) Two typical examples of the 
estimated PRCs, (b) The dependence of the PRC shape on 
a. (c) The distributions of the estimated shape parameters 
(upper) and the intrinsic frequencies of neuronal oscillators 
(lower) in different cortical layers. The frequency was tuned 
in experiments by varying the amplitude of a DC injected 
current. 



We first derive the interaction functions Tj (if>) [|| from 
a set of differential equations of the phase oscillators with 
a common frequency u). These oscillators have heteroge- 
neous PRCs, Zj(0), and are globally coupled with each 
other. The phase of the jth oscillator 8j obeys the evo- 
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lution equation, 
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The input to the jth oscillator from the fcth is given as 
eN~ x J2 n a (t ~ *fc)> w here e is a weak coupling constant, 
N the number of oscillators, a(t) a causal coupling func- 
tion, and the nth firing time of the kth oscillator fZ de- 
fined as Ok(tfc) — 27m. 

The mutual interactions shift the frequency of the 
mean phase of these oscillators by eCl from the natural 
frequency u. We define the relative phase ipj = 9j — $, 
where the phase $ = (u> + ef2)t. The relative phase ipj 
changes slowly compared with 9j and will hardly change 
during the oscillation period 2n/{u + eft). Therefore, we 
can average Eq. (TT]) over one period keeping ipj constant. 
Using (id + e£l)t k = 2nn — tpk, we can describe the dy- 
namics of the relative phase ipj as 
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where f3(ipk + <&) is the sum of the contributions of 
past inputs from the fcth oscillator ^ n ct{{ipk + $ — 
2-nn)l{io + eft)). In the limit of N — > oo, we can ap- 
ply the mean-field approximation, iV -1 J2k Pii'k + <&) ~ 
Ju* j3(ip + $>)P(ip)dip, where P(ip) is the distribution 
function of the relative phase ip. The functions Z, P and 
(3 are 27r-periodic, so we can expand them into Fourier 
series. The nth Fourier coefficient Tn ' l e lX ™ > of a function 
f(ip) is defined as r n f) e ixi ^ = (27T)- 1 f(x)e' mx dx. 
Then, Eq. ([2]) becomes, 

e dt jW 
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the phase distribution function of the 'synchronized' os- 
cillators, P s , we can represent 0™ as 
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where D refers to the indices of the 'synchronized' os- 
cillators and ipj is a real- valued solution to T(ipj) = 
satisfying T'(ipj) < 0. In general, the equation has more 
than one solution. However, Eq. Q can uniquely be de- 
fined in the limit of weak noise since the noise excludes 
solutions other than the most stable one. 

The contribution of the 'desynchronized' oscillators 
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where D is the set of the indices of the 'desynchronized' 
oscillators and Pds is their phase distribution function. 
We can calculate the above integral using the residue 
theory as 
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be- 



ing an imaginary solution to T(ip^) — satisfying 
Im ipj* < and the weighted average (f(xk);g(xk)) k 

defined as (f(x k );g(x k )) k = J2k g{xk)f(xk)/ J2k 9( x k)- 
We find that the formal expressions of the contributions 
of the 'desynchronized' and 'synchronized' populations 
are identical. Thus, we finally obtain the following self- 
consistent equation: 



where U)j = r^rf*, Kf = 2r ( n j) r { n ] and A™ = - 

\n . The order parameters R n = 2-Kr n F ' , fl, and 
in Eq. describe the dynamics of neural oscillators. 

Below, we focus on non-trivial solutions (R n ^ 0) to 
Eq. ([3]) other than the trivial one R n = 0. Each oscillator 
exhibits two dynamical modes, 'synchronized' or 'desyn- 
chronized', according to the shape of the PRC or the 
interaction function Tj(ip). In the 'synchronized' mode, 
the oscillator is trapped at a stable fixed point of Eq. ([3]) . 
In the 'desynchronized' mode, Eq. © has no stable fixed 
point, so the oscillator cannot be locked at any relative 
phase and drifts at a period of Tj = (T 1 ^ \Tj (ip) \~ 1 dip. 
To derive the order-parameter equations, we introduce 



complex order parameters R n e 



;a< p > 



and divide them into 



the contributions of the 'synchronized' population, 0™, 
and those of the 'desynchronized' population, 0^ s . Given 
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This equation means that the nth complex order param- 
eters i?™e lA " ' should be identical with the nth circular 
moment of the complex solutions to Tj(ip) = 0. 

To obtain the explicit formula for the fixed points, we 
hereafter truncate T(ip) up to the first Fourier mode of 
Eq. ©. Then, the complex solutions ip* are given by 
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(Cjj — Cfy/KjR 1 . If a single variable a pa- 
rameterizes the heterogeneity of the PRC shapes, we can 
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explicitly describe the phase distribution functions as 



P B (i/>) = g(a{xl>)) 
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where g(a) is the distribution of the parameter values 
over the oscillator population. As mentioned previously, 
T(ip,a) — has no real solution in the parameter range 
-Dds- It is noted that merely finding the equilibrium val- 
ues of the order parameters does not require the explicit 
expressions of P s ^ s (ip) m the present analysis. 

To show the validity of our theoretical treatment and 
to get a novel insight into the dynamics of heterogeneous 
oscillator networks, we now apply it to coupled oscillators 
having the following heterogeneous PRCs: 
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an) + cos an. 
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The value of the shape parameter a is distributed uni- 
formly in the range (a nl - m ^ a ^ a max )- A large or 
a small value of a corresponds to type II- or type I- 
like PRC, respectively. The coupling function ot(ip) is 
an exponential function with a decay constant of r: 
a(t) = T~ 1 Q(t)e~ t / T , where 0(x) is Heaviside function: 
0(a;) = 1 if x > or if x < 0. Neurons synchronize 
easier with type II-like PRCs than with type I- like PRCs. 
Using Eq. (fTTj) . we can rewrite Eq. ((3|) as 
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W(a) 



W(a) — cos (ipj — an + Arctan(rw)) 
cosa7r — 2nfl/uj 



(12) 



where K\ = w/(2vr^l + (tw) 2 ). 

Using these results, we can study different dynamical 
states of the oscillator network. Figure Ufa) summarizes 
the phase diagram in the (a mm , a max ) half plane. When 
Eq. (fT3)) has a non-trivial solution, neurons are either 
partially synchronized (PaS) or perfectly synchronized 
(PfS). The border between the two states by can be de- 
termined by a critical value a c , which is a solution to 
|W(a c )| = 1. Since neurons with a > a c are synchro- 
nized and those with a < a c are desynchronized, the PfS 
state requires a m - m > a c . Otherwise, neurons are only 
partially synchronized. If a m [ n is sufficiently large, all 
neurons may be type II-like. We note that even in such 
a case, a strong heterogeneity (i.e., a sufficiently large 
flmax — Omin) may disable the perfect synchronization of 
oscillators. 

The entire population of oscillators is perfectly desyn- 
chronized (PfD) if the self-consistent equation, 



R 1 
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(13) 



has no non-trivial (R 1 ^ 0) solutions in the allowed range 
of et-value. Thus, the border between the PaS and PfD 



states is determined from the condition that a solution 
with R 1 ^ exists. At a m i n = a max , the self-consistent 
equation has a stable fixed point if a > 7r~ 1 Arctan(ru;) 
and the system exhibits the PfS state (point Q). The PfS- 
PaS and PaS-PfD borders should merge at this point. 
Figure [2](c)-(e) displays the raster plots of the neural 
oscillators in the PfS, PaS, and PfD states designated 
in Figure [DJa), respectively. In Figure [^c) or (e), the 
population comprises only type II-like or type I-like neu- 
rons showing perfect synchronization or perfect desyn- 
chronization, respectively. In Figure [21(d) , only sub- 
population of strongly type II-like neurons with a > a c 
are synchronized. 
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FIG. 2: Three different states in our coupled oscillator model 
with the PRCs represented by Eq. . The parameters were 
set as t = 0.005, id = 40n, and e = 0.01. (a) Phase diagram 
of this model, where the shape parameter a is distributed uni- 
formly in o m in ^ = a m ax- The abbreviations mean perfectly 
synchronized (PfS), partially synchronized (PaS), and per- 
fectly desynchronized (PfD) states, (b) The order parameters 
R 1 and a m i n — a c are shown along the line a max — a m i n = 0.1 
(dashed line in (a)). The PfD state corresponds to R 1 = 0, 
while the PaS or the PfD state is defined with R 1 7^ and 
with Omin — a c < or a m in — a c > 0, respectively, (c)— (e) 
Raster plots of the relative phases ipj of 100 neural oscillators 
in the PfS (c), PaS (d), and PfD (e) states. 



The necessary condition for getting a non-trivial solu- 
tion to Eq. (fl3"|) is that the imaginary part of its right 
hand should vanish. While the original Kuramoto model 
always satisfies this condition due to its symmetry, this 
is not the case in the present example. Therefore, this 
model and the Kuramoto model exhibit qualitatively dif- 
ferent transitions from the PfD to the PaS state. In this 
model, the order parameter R 1 , which vanishes in the 
PfD state, jumps discontinuously to a non-zero value at 
the transition point (Fig. HJb)). In contrast, the tran- 
sition is continuous in the Kuramoto model. Therefore, 
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FIG. 3: The layer dependences of the PRCs of cortical neu- 
rons in the 7- frequency range and the order parameter R 1 . 
(a) The PRCs recorded in pj| were fitted by Eq. (JTTJ. PRCs 
were superimposed for 10 layer-II/III neurons (broken black 
lines) and 25 layer- V neurons (solid gray lines), (b) The or- 
der parameter R 1 depends differentially on the frequency us 
of neuronal firing for layer-II/III (black) and layer- V (gray) 
neurons. Here, r = 0.005. Eq. (|13p in general has multiple 
solutions when the summation is taken over a finite number 
of oscillators. Here, we plotted the largest R , or plotted 
nothing if R 1 < 0.1. Dashed line indicates a discontinuous 
transition of the network state. Other discontinuous points 
of the curves appeared due to the finite-size effect. 



the network state switches very sharply in the present 
model. Our quantitative results reveal that the differ- 
ent types of the heterogeneity can result in qualitatively 
different phase transition-like behaviors of oscillator net- 
works. 

We apply our theory to the data recorded from excita- 
tory neurons in cortical layers II/III and V [15!]. The 



PRC shapes exhibited a remarkable layer dependence 
when the firing frequency w is in a range of 20-45 Hz, 
which is within the 7- frequency range (Fig. EJa)). To 
demonstrate the effect of this layer dependence on the 
population dynamics of cortical neurons, we analyze cou- 
pled systems of layer-II/III or layer- V neurons separately. 
The synchronizing property in general depends on the 
frequency us as well as the decay time constant of excita- 
tory synapses. As us is increased, layer-V neurons display 
an abrupt transition from the PaS to PfD state, as repre- 
sented by a discontinuous jump in R 1 (dashed line). This 
transition can also be found in the simpler model shown 
previously in Fig. [2] with a uniform distribution of a. 

These results may have considerable implications for 
exploring computational functions of local cortical cir- 
cuits. Furthermore, in most real-world oscillator systems, 
the phase response curves are highly heterogeneous, so 
our theoretical method can be naturally applied to these 
systems. Our results provide a way to make quantitative 
predictions about the dynamics of general heterogeneous- 
oscillator networks, indicating the existence of a clear-cut 
border between the PaS and PfD states with a discontin- 
uous jump in the order parameter. 

We thank Y. Kuramoto, T. Mizuguchi, H. Nakao, 
G. Bi, H. Cateau and N. Masuda for fruitful discus- 
sions and valuable comments. This work has been sup- 
ported in part by Grant-in-Aid for Scientific Research 
in Priority Areas (17022036) and for Young Scientists 
(B) (50384722) from the Japanese Ministry of Education, 
Culture, Sports, Science and Technology. 



[1] P. Hadley, M.R. Beasley, and K. Wiesenfeld, Phys. Rev. 

B 38, 8712 (1988); S. Watanabe and S.H. Strogatz, Phys. 

Rev. Lett. 70, 2391 (1993). 
[2] S.S. Wang and H.G. Winful, Appl. Phys. Lett. 52, 1774 

(1998). 

[3] D. Plenz and S.T. Kitai, J. Neurophysiol. 76, 4180 
(1996); S.F. Farmer, J. Physiol. 509, 3 (1998); P. Fries, 
J.H. Schroder, P.R. Roelfsema, W. Singer, and A.K. En- 
gel, J. Neurosci. 22, 3739 (2002); C.S. Herrmann, M.H. 
Munk, and A.K. Engel, Trends Cogn. Sci. 8, 347 (2004); 
D. Lee, J. Neurosci. 24, 4453 (2004). 

[4] A.T. Winfree, J. Theor. Biol. 16, 15 (1967). 

[5] A.T. Winfree, The Geometry of Biological Time 
(Springer, New York, 1980). 

[6] Y. Kuramoto, Chemical Oscillations, Waves, and Turbu- 
lence (Springer- Verlag, Berlin, 1984). 

[7] A.D. Reyes and E.E. Fetz, J. Neurophysiol. 69, 1661 
(1993); A.D. Reyes and E.E. Fetz, J. Neurophysiol. 69, 
1673 (1993). 

[8] D. Hansel, G. Mato, and C. Meunier, Neural Comput. 7, 
307 (1995). 

[9] B. Ermentrout, Neural Comput. 8, 979 (1996). 



[10] R.F. Galan, G.B. Ermentrout, and N.N. Urban, Phys. 

Rev. Lett. 94, 158101 (2005). 
[11] A.J. Preyer and R.J. Butera, Phys. Rev. Lett. 95, 138103 

(2005) . 

[12] M. Lengyel, J. Kwag, O. Paulsen, and P. Dayan, Nat. 

Neurosci. 8, 1667 (2005). 
[13] T.I. Netoff, M.I. Banks, A.D. Dorval, CD. Acker, J.S. 

Haas, N. Kopell, and J. A. White, J. Neurophysiol. 93, 

1197 (2005). 

[14] T. Tateno, and H.P. Robinson, Biophys. J. 92, 683 
(2007). 

[15] Y. Tsubo, M. Takada, A.D. Reyes, and T. Fukai, Eur. J. 

Neurosci. 25, 3429 (2007). 
[16] H. Sakaguchi and Y. Kuramoto, Prog. Theor. Phys. 76, 

576 (1986). 

[17] T. Shimokawa, S. Shinomoto, Phys. Rev. E. 73, 066221 

(2006) . 

[18] H. Daido, Phys. Rev. Lett. 87, 048101 (2001). 

[19] J. Teramae and D. Tanaka, Phys. Rev. Lett. 93, 204103 

(2004); J. Teramae and D. Tanaka, Prog. Theor. Phys. 

Supp. 161, 360 (2006). 



